#########################################################################
# Replication for JPR Resolving Bargaining Problems in Civil Conflict   #
# Joo 2024 						   									                              #
# Models run on Dell Latitude 7420										                  #
# Windows 10 Enterprise (V.22H2)                                        #
# Processor: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz   1.80 GHz  # 
# RAM: 16.0 GB (15.7 GB usable)    										                  #
# R version 4.3.2 (2023-10-31 ucrt) -- "Eye Holes"                      #
#########################################################################

rm(list=ls())
library(MASS)
library(mvtnorm)
library(ggplot2)
library(foreign)

wd <- setwd("Your Directory/JPR_Joo_replication2024")

negdata <- read.dta("nego_maindata.dta")

set.seed(20240129)

coef <- c(-0.2277705,
          0.1462724,
          0.7675535,
          -0.1892973,
          -0.0578004,
          0.2862604,
          -0.2117465,
          -0.4033382,
          0.041031,
          -0.0011928,
          -0.3245317)

vcov <- read.csv("vcov.csv")

vcov <- as.matrix(vcov)

is <- 0 
com <- 0 
tc <- 0  
pw <- 0
med <- 0  
rsup <- 1
gsup <- 1 
war <- 0
bd <- mean(negdata$lbdbest, na.rm =T) 
par <- 0
t1 <- mean(negdata$T1, na.rm=T)
t2 <- mean(negdata$T2, na.rm=T)
t3 <- mean(negdata$T3, na.rm=T)

poly <- seq(0,1, by=.01)
set.seed(20240129)
B = 1000


#############
# Figure 1a # 
#############


pr.bs <- matrix(NA,B,2)
SD <- sd(negdata$polyarchy, na.rm=T)
edemo <- mean(negdata$polyarchy, na.rm=T)


set.seed(20230905)

for(i in 1:B){ 
  newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
  
  ### moderate = 1 ###
  
  # Edem 1 SD below mean    
  bs.sum1 <- newcoef[1]*1 + newcoef[2]*(edemo-SD) + newcoef[3]*1*(edemo-SD) + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
    newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
  
  # Edem 1 SD above mean 
  bs.sum2 <- newcoef[1]*1 + newcoef[2]*(edemo+SD) + newcoef[3]*1*(edemo+SD) + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
    newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
  
  pr.bs[i,1] <- pnorm(bs.sum2) - pnorm(bs.sum1)
  
  
  ### moderate = 0 ###
  
  # Edem 1 SD below mean 
  bs.sum3 <- newcoef[1]*0 + newcoef[2]*(edemo-SD) + newcoef[3]*0*(edemo-SD) + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
    newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
  
  
  # Edem 1 SD above mean 
  bs.sum4 <- newcoef[1]*0 + newcoef[2]*(edemo+SD) + newcoef[3]*0*(edemo+SD) + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
    newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
  
  pr.bs[i,2] <- pnorm(bs.sum4) -  pnorm(bs.sum3) 
  
}

mod1 <- data.frame(results=pr.bs[,1], model="mod1")
m1 <- as.numeric(mean(mod1$results))
u1 <- as.numeric(quantile(mod1$results, p=0.975))
l1 <- as.numeric(quantile(mod1$results, 0.025))
u12 <- as.numeric(quantile(mod1$results, p=0.995))
l12 <- as.numeric(quantile(mod1$results, 0.005))


mod0 <- data.frame(results=pr.bs[,2], model="mod0")
m0 <- as.numeric(mean(mod0$results))
u0 <- as.numeric(quantile(mod0$results, 0.975))
l0 <- as.numeric(quantile(mod0$results, 0.025))
u02 <- as.numeric(quantile(mod0$results, 0.995))
l02 <- as.numeric(quantile(mod0$results, 0.005))

CIdata <- data.frame(model=c("Maximalist", "Moderate"), mean=c(m0, m1), 
                     lower=c(l0, l1), upper=c(u0,u1), lower2 = c(l02, l12), upper2=c(u02, u12))

jpeg(paste(wd, "/Figures/Figure1a.jpeg", sep=""), width=6,height=6,units="in", res=400)

par(mar=c(3,2,1,1))

ggplot() + 
  geom_errorbar(data=CIdata, mapping=aes(x=model, ymin=lower, ymax=upper),  width=0.2, linewidth=1, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=22),
                     axis.title.x = element_text(size=18),
                     axis.text.x = element_text(size=18),
                     axis.text.y = element_text(size=14))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2.5,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  linewidth=1) + 
  xlab("Rebel Goal") + ylab("Change in Pr(Negotiation=1)") + ylim(-0.12, 0.21)

dev.off()

#############
# Figure 1b #
#############


out2 <- matrix(NA, ncol=3, nrow=length(poly))

for(i in 1:length(poly)){ 
  edem <- poly[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*0 + newcoef[2]*edem + newcoef[3]*0*edem + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
      newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
    
    bs.sum2 <- newcoef[1]*1 + newcoef[2]*edem + newcoef[3]*1*edem + newcoef[4]*is + newcoef[5]*com +newcoef[6]*tc + 
      newcoef[7]*pw + newcoef[8]*t1 + newcoef[9]*t2 + newcoef[10]*t3 + newcoef[11]
    
    pr.bs[b,] <- plogis(bs.sum2) - plogis(bs.sum) 
  }
  out2[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

jpeg(paste(wd, "/Figures/Figure1b.jpeg", sep=""), width=6,height=6,units="in", res=400)

par(mar=c(5,5,1,1))
plot(poly, out2[,1], type="n", xlab="", ylab="", ylim=c(-0.1, 0.23), cex.axis=1)
lines(lowess(poly, out2[,1]), lty=1, lwd=2)
lines(lowess(poly, out2[,2]), lty=2, lwd=2)
lines(lowess(poly, out2[,3]), lty=2, lwd=2)
abline(h=0, lty=2, col='red', lwd=2)
title(xlab="Polyarchy", line = 3, cex.lab = 1.7)            
title(ylab="Change in Pr(Negotiation=1)", line = 3, cex.lab=1.7)  

dev.off()


# End of file